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AN  IMPROVED  LIMITER  FOR 
MULTIDIMENSIONAL  FLUX-CORRECTED  TRANSPORT 


1.  INTRODUCTION 

Flux-corrected  transport,  or  FCT  (Boris  Si  Book  1973;  Book,  Boris  Si  Hain  1975;  Boris 
Si  Book  1976a,  1976b;  Zalesak  1979;  DeVore  1989, 1991;  Boris  et  al.  1993),  was  developed, 
and  has  been  used  extensively  by  many  investigators,  to  accurately  solve  the  conservation 
equations  of  Eulerian  hydrodynamics  without  violating  the  positivity  of  mass  and  energy, 
particularly  in  the  vicinity  of  shock  waves  and  other  discontinuities.  This  is  achieved  by 
adding  to  the  equations  a  strong  numerical  diffusion,  which  guarantees  the  positivity  of 
the  solution,  followed  by  a  compensating  antidiffusion,  which  reduces  the  numerical  error. 
The  crux  of  the  FCT  method  lies  in  limiting  (‘correcting’)  these  antidiffusive  fluxes  before 
they  are  applied,  so  that  no  unphysical  extrema  are  created  in  the  solution.  The  goal  of 
the  flux-correction  procedure  is  to  provide  as  accurate  a  solution  to  the  original  equation 
as  is  consistent  with  maintaining  positivity  and  monotonicity  everywhere. 

The  original  limiter  of  Boris  and  Book  (1973)  applies  several  tests  to  the  sign  and 
magnitude  of  the  antidiffusive  flux  at  each  cell  boundary,  taking  into  account  the  profile 
of  the  transported  quantity  (e.g.,  the  mass  or  energy  density)  in  the  neighborhood  of  that 
boundary.  In  one  spatial  dimension,  these  tests  can  be  condensed  into  a  single,  relatively 
simple  formula  for  the  corrected  antidiffusive  fluxes,  which  then  are  guaranteed  to  preserve 
the  positivity  and  monotonicity  of  the  transported  variable.  For  certain  multidimensional 
hydrodynamics  applications,  this  limiter  and  its  underlying  integration  scheme  can  be 
employed  in  a  serial  fashion,  to  each  of  the  coordinate  directions  in  turn,  using  Strang 
operator  splitting  (Oran  Si  Boris  1987). 

Certain  classes  of  problems,  however,  call  for  a  fully  multidimensional  approach.  These 
include  incompressible  or  nearly  incompressible  flows,  flows  with  a  high  degree  of  symmetry, 
and  magnetohydrodynamic  (MHD)  systems.  Zalesak  (1979)  analyzed  Boris  and  Book’s 
one-dimensional  limiter,  identified  the  tests  inherent  in  it,  and  showed  why  a  reformulation 
was  needed  to  extend  the  FCT  approach  to  multidimensional  situations.  He  provided  such 
a  formulation,  showed  that  his  new  limiter  could  be  made  equivalent  to  the  original,  and 
then  illustrated  it  with  several  one-  and  two-dimensional  examples. 

Zalesak’s  multidimensional  limiter  has  been  used  by  the  author  in  both  hydrodynamic 
and  MHD  (DeVore  1989;  DeVore  1991)  FCT  algorithms  in  recent  years,  with  applications 
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to  shock  and  detonation  waves  (Oran  k  DeVore  1994),  shear  flows  and  jets  (Grinstein  1995; 
Grinstein  k  DeVore  1996),  and  arcades  of  plasma  and  magnetic  field  in  the  solar  atmo¬ 
sphere  (Karpen,  Antiochos  k  DeVore  1990,  1991,  1994,  1996;  Karpen  et  al.  1998;  Antio- 
chos,  DeVore  k  Klimchuk  1999;  Antiochos  k  DeVore  1999a,  1999b),  among  others.  Those 
experiences  have  shown  that  the  limiter  is  positivity-  but  not  monotonicity-preserving, 
yielding  solutions  with  numerical  ripples  of  significant  amplitude  in  many  cases.  The  pur¬ 
pose  of  this  paper  is  to  show  why  this  occurs  and  to  propose  a  modification  of  Zalesak’s 
limiter  which  contributes  greatly  towards  preserving  monotone  profiles.  The  recommended 
modification  is  the  addition  of  a  prelimiting  step  based  on  Boris  and  Book’s  original  limiter, 
which  is  both  positive  and  monotone.  Zalesak  in  fact  suggested  this  approach  in  passing 
in  his  earlier  work,  but  did  not  explore  its  ramifications  and  did  not  promote  its  regular 
use.  We  do  both  in  the  remainder  of  this  paper. 

2.  GENERALIZED  CONTINUITY  EQUATIONS  AND  FCT 

The  conservation  equations  of  Eulerian  hydrodynamics  take  the  general  form 

^  +  V  •  (pv)  =  S,  (1) 

where  p  is  a  fluid  variable  (mass,  momentum,  or  energy  density)  being  time-advanced,  v 
is  the  fluid  velocity,  and  S  is  a  source  term.  In  a  finite-volume,  Eulerian  representation 
of  Eq.  (1),  the  mass,  momentum,  or  energy  enclosed  in  each  computational  cell  evolves 
due  to  sources  and  due  to  the  convective  fluxes  through  the  faces  of  that  cell.  Each  such 
flux  defined  on  an  interior  face  of  the  domain  increments  the  amount  of  mass  in  the  cell 
which  the  flux  is  entering,  and  decrements  it  by  an  equal  amount  in  the  cell  which  it  is 
exiting.  This  leaves  the  total  mass  of  the  system  unchanged,  i.e.,  the  integration  scheme 
is  conservative. 

Any  discrete  representation  of  the  convective  fluxes  in  Eq.  (1)  introduces  numerical 
errors  into  the  solution  for  p.  Indeed,  a  central  difference  representation  of  the  fluxes  - 
wherein  the  face  value  of  p  is  just  the  average  of  the  cell  values  ahead  of  and  behind  the 
face  -  produces  an  absolutely  unstable  solution  for  any  finite  velocity  v.  These  errors  can 
be  reduced  by  introducing  additional  numerical  terms  into  Eq.  (1),  which  also  are  defined 
in  terms  of  fluxes  so  that  the  integration  scheme  remains  conservative.  In  particular, 
the  addition  of  a  smoothing  diffusion  term  of  sufficient  amplitude  stabilizes  the  convective 
terms  and  also  imposes  one  other  important  physical  constraint  on  the  solution:  positivity. 
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This  ensures  that  the  evolution  of  a  positive-definite  initial  profile  p°  does  not  lead  to 
unphysically  negative  densities  as  a  consequence  of  convection  alone. 

The  result  of  this  integration  procedure  is  a  provisional  solution  pl  which  is  positive 
except  where  the  sources  S  are  negative  and  sufficiently  large  to  produce  a  reversal  in 
sign,  and  is  monotone  except  where  either  convection  or  sources,  or  both  together,  act  to 
produce  extrema. 

The  strong  numerical  diffusion  required  to  guarantee  the  positivity  of  the  solution  (f 
is  too  large  to  be  satisfactory  for  most  hydrodynamics  applications.  In  regions  where  the 
density  is  slowly  varying,  it  is  advantageous  to  apply  antidiffusive  fluxes  which  compensate 
for  the  previously  applied  diffusion  and  yield  a  high-order  accurate  solution  ph  to  Eq. 
(1).  Unfortunately,  near  discontinuities  this  procedure  is  perilous,  as  the  attempt  to  fit 
a  high-order  solution  can  produce  a  highly  oscillatory  solution.  Both  the  positivity  and 
monotonicity  properties  of  the  provisional  solution  can  be  lost  in  this  effort. 

With  flux-corrected  transport  (FCT),  Boris  and  Book  (1973)  resolved  the  conflicting 
demands  for  a  positive,  monotone  solution  such  as  //  on  the  one  hand,  and  for  an  accurate, 
nondiffusive  solution  such  as  ph  on  the  other.  Their  approach  is  to  limit  the  antidiffusive 
fluxes  which  change  pi  into  ph  so  that  no  new  extrema  are  created ,  nor  are  any  existing 
extrema  enhanced  by  the  application  of  the  corrected  fluxes.  The  rule  which  they  imple¬ 
mented  to  achieve  this  objective  will  be  considered  momentarily.  The  FCT  prescription 
for  solving  generalized  continuity  equations  such  as  Eq.  (1)  can  be  summarized  as  follows: 

1.  Convect  p. 

2.  Numerically  diffuse  p  sufficiently  to  ensure  that  the  result  is  positive  everywhere,  if 
the  initial  profile  p°  is  positive. 

3.  Add  sources  S  to  obtain  the  low-order,  provisional  solution  pl. 

4.  Compute  numerical  antidiffusive  fluxes  which  would  convert  pe  to  a  high-order  accu¬ 
rate  solution  ph. 

5.  Limit  these  antidiffusive  fluxes  so  that  no  new  extrema  will  be  created  and  no  existing 
extrema  will  be  enhanced. 

6.  Apply  the  corrected  fluxes  to  p£  to  obtain  the  desired  solution  pn . 

The  combination  of  the  positive,  monotone  solution  from  steps  1  and  2  with  the  restrictive 
action  of  the  flux  limiter  in  step  5  ensures  that  the  result  pn  also  will  be  positive  and 
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Fig.  1.  The  action  of  Boris  and  Book’s  flux  limiter  for  one-dimensional  FCT  is  indicated  for  the  eight 
possible  configurations  of  four  grid  points  (•)  surrounding  an  antidiffusive  flux  (—>)  here  shown  directed 
to  the  right.  The  dashed  lines  in  the  leftmost  panel  indicate  the  maximum  excursion  allowed  the  two 
middle  points  when  the  corrected  flux  is  applied,  so  that  no  new  extrema  are  created.  In  the  middle  and 
rightmost  panels,  the  flux  is  cancelled  to  avoid  enhancing  existing  extrema.  On  the  bottom  row,  the  flux 
is  reversed  in  direction  so  that  it  acts  to  steepen  the  profile. 

monotone  except  where  the  physics  of  convection  and  sources  together  produce  extrema. 
Spurious  numerical  extrema  are  avoided. 


3.  BORIS  AND  BOOK’S  LIMITER 


In  their  original  paper  on  FCT,  Boris  and  Book  (1973)  provided  a  rule  for  limiting  the 
antidiffusive  fluxes  when  solving  generalized  continuity  equations  in  one  spatial  dimension. 
Given  the  provisional  profile  pf,  the  cell  volumes  Vj ,  and  the  raw  antidiffusive  fluxes 
defined  on  the  cell  faces,  where  i  is  a  cell  index,  the  rule  is 


F?+ 1/2  =  cri+ 1/2  max  0,  min  (  F?+1/ 2  ,  <?i+ 1/2  Vi  Apf_1/2,  <r*+i/2  Vi+i  ^Pi+3/2)  _ 


(2) 
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Here  F?+1/2  is  the  corrected  (positivity-  and  monotonicity-preserving)  antidiffusive  flux, 
Ap£  is  a  density  difference, 

aP?+i/2  =Pi+i  -Pi>  (3) 

and  a  is  its  sign, 

^i+i/2  =  signApf+1/2,  (4) 


with  |cr|  =  1. 

This  flux  limiter  has  the  following  consequences: 

•  If  the  pf:  profile  is  not  monotone  between  cells  i- 1  and  *+ 2,  which  means  that  at  last 
one  of  pf  and  pei+1  is  a  local  extremum  and  A p£  changes  sign,  then  the  flux  is  cancelled. 

•  If  the  pl  profile  is  monotone  over  that  range,  then  the  amplitude  of  the  flux  is  limited 
to  the  smallest  of  its  original  magnitude  and  the  two  mass  changes  sufficient  to  flatten 
the  profile  ahead  of  and  behind  the  interface. 

•  If  the  raw  antidiffusive  flux  is  directed  down  the  gradient  of  pe,  so  that  it  smooths 
rather  than  steepens  the  profile,  then  the  sign  of  the  flux  is  reversed  (it  is  given  the 
sign  of  Ap^)  and  the  amplitude  is  limited  according  to  the  above  criteria. 

These  consequences  are  displayed  schematically  in  Fig.  1. 

A  consideration  of  the  flux  limiter  (2)  and  its  actions  in  the  various  scenarios  sketched 
in  Fig.  1  shows  that  its  stated  goal  is  met:  no  new  extrema  are  created,  nor  are  any 
existing  extrema  enhanced  as  a  result  of  the  antidiffusion  stage  of  the  algorithm.  The 
former  is  ensured  by  the  amplitude  limitation  inherent  in  the  flux-correction  formula;  the 
latter  by  the  complete  cancellation  where  the  density  difference  reverses  sign. 

4.  ZALESAK’S  LIMITER 

In  his  analysis  of  the  flux  limiting  formula  (2)  to  derive  an  extension  of  it  to  two  or 
more  spatial  dimensions,  Zalesak  (1979)  noted  the  remarkable  fact  that  each  flux  is  limited 
independently;  there  is  no  consideration  of  the  effect  of  multiple  fluxes  acting  in  concert. 
Upon  reflection,  this  is  to  be  expected  in  one  dimension  because  both  fluxes  enter  or  leave 
a  cell  only  if  that  cell  is  a  local  maximum  or  minimum,  respectively.  In  either  instance, 
both  fluxes  must  be  cancelled  to  avoid  accentuating  the  already  existing  extremum.  This 
cancellation  is  already  built  into  the  flux  limiter,  however,  so  knowledge  about  the  sign 
and  magnitude  of  neighboring  fluxes  is  not  required  in  (2). 
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The  situation  is  different  in  two  or  more  dimensions,  where  multiple  fluxes  may  enter 
or  leave  a  cell  without  that  cell  being  a  local  extremum.  As  a  result,  the  flux  limiter 
must  take  into  account  the  consequences  of  fluxes  acting  in  concert.  Zalesak  reformulated 
the  one- dimensional  flux  limiter  in  a  way  that  preserves  the  desirable  properties  of  Boris 
and  Book’s  limiter,  but  also  generalizes  immediately  to  multidimensional  problems.  His 
procedure  is  as  follows: 

A.  Establish  allowed  extrema  pmin  and  pmax  in  each  cell;  for  consistency  with  Boris  and 
Book,  set 

pf“  =  min  (pi-i>Pi> Pi+i)  >  ^ 

P?AX  =  max  {pi-i,  pi  Pi+i)  ■ 

B.  Reverse  the  sign  of  any  antidiffusive  flux  that  is  directed  down  the  gradient  of  pe,  viz., 


Fi+ 1/2  —  °*+ 1/2 


T?a 

*i+l/2 


(6) 


C.  Compute  the  total  incoming  and  outgoing  antidiffusive  fluxes  Fm  and  Fout  in  each 
cell, 

Ef  =  max  (i?:i/2,0)  -  min  (f£1/2,  OJ  , 

Ffut  =  max  (F?;i/2,0)  -  min  (f?1  i/2,o)  . 

D.  Determine  the  fractions  of  the  incoming  and  outgoing  fluxes  that  can  be  applied  to 
each  cell, 

ft  =  (pr*  -  $  Vi/ 

St1  =  (Pi  -  Pfto)  V,/FtM. 

E.  Limit  each  antidiffusive  flux  so  that  it  creates  neither  an  undershoot  in  the  cell  it  is 
leaving  nor  an  overshoot  in  the  cell  it  is  entering, 


*tn/2  =  *7+1/2  >< 


min  (/fut,  /l+n  l)  >  if  ^T+i/2  ^  °5 
min  (fln ,  ,  1 ) ,  otherwise. 


(9) 


Zalesak’s  step  B  as  given  in  his  paper  is  different  from  the  one  above,  and  his  is  not 
fully  consistent  with  the  original  limiter.  The  discrepancy  is  that  he  used  the  sign  of  the 
antidiffusive  flux  Fa,  rather  than  the  sign  of  the  density  difference  A pe,  for  a  in  Eq.  (2). 
As  he  noted,  this  adjustment  is  relatively  rarely  applied  in  any  case  -  the  majority  of 
raw  antidiffusive  fluxes  act  to  steepen  the  gradient  -  but  it  does  serve  to  suppress  the 
generation  of  numerical  ripples.  Also,  Zalesak  pointed  out  and  showed  through  several 
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Fig.  2.  A  quasi-one-dimensional  discontinuity  propagating  in  2D  is  shown  as  a  surface  plot.  Each  junction 
of  lines  represents  a  grid  point  (cell  center)  where  the  density  p *  is  defined.  Transverse  to  the  principal 
discontinuity,  which  is  oriented  left  to  right,  is  a  ripple  in  the  solution.  Zalesak’s  flux  limiter  in  2D  will 
allow  this  ripple  to  be  created  and  amplified,  so  long  as  the  extrema  along  the  ripple  lie  within  those  of 
the  principal  discontinuity  transverse  to  it. 

examples  how  the  limiter  can  be  made  more  flexible  by  generalizing  the  formulae  of  step 
A  for  the  allowed  extrema. 

The  flux  limiter  as  reformulated  by  Zalesak  has  the  following  consequences: 

•  It  can  be  made  identical  to  Boris  and  Book’s  limiter  in  one  spatial  dimension. 

•  It  generalizes  readily  to  multiple  dimensions. 

•  It  enforces  positivity  if  pl  is  positive  and  if  the  allowed  extrema  pmm  and  pmax  are 
adequately  constrained. 

•  It  prevents  the  creation  of  new  extrema  and  the  enhancement  of  existing  ones. 

•  However,  the  creation  of  new  and  the  enhancement  of  existing  numerical  ripples  in  the 
solution  pn  is  allowed  in  multidimensional  problems,  i.e.,  the  limiter  is  not  monotone. 
To  see  that  these  last  two  statements  are  not  in  conflict,  consider  the  quasi-one- 

dimensional  discontinuity  shown  schematically  in  Fig.  2.  Imposed  on  the  steep  portion 
of  the  profile  is  a  transverse  ripple  which  the  high-order  convection  algorithm  can  create 
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initially  or  amplify  from  its  current  amplitude,  above  and  beyond  what  the  convection  and 
sources  may  produce  physically.  Zalesak’s  limiter  does  not  prevent  this  from  happening 
so  long  as  the  excursions  associated  with  the  ripple  do  not  surpass  the  local  extrema. 
These  extrema  are  located  ahead  of  and  behind  the  ripple,  along  the  direction  of  principal 
variation  in  the  profile.  The  ripple  shown  may  grow  until  its  maximum  or  minimum  value 
reaches  those  of  the  neighboring  extrema. 


5.  A  MONOTONE  MULTIDIMENSIONAL  LIMITER 

Clearly,  it  is  desirable  that  the  transverse  ripple  in  Fig.  2  not  be  amplified  during 
the  antidiffusion  stage  of  the  FCT  solution,  nor  should  it  be  created  in  the  first  place, 
unless  the  convection  and  sources  conspire  to  produce  it.  This  suggests  that  the  limiter 
be  extended  to  prevent  the  creation  of  new  and  the  enhancement  of  existing  directional 
extrema,  i.e.,  extrema  with  respect  to  each  coordinate  direction  considered  independently. 
The  solution  adopted  here  is  to  apply  Boris  and  Book  s  limiter  along  each  coordinate 
direction  separately,  as  a  prelimiting  step  in  the  FCT  procedure.  Within  the  context  of 
Zalesak’s  limiter,  this  prelimiting  step  replaces  and  augments  the  sign-reversal  step  B.  It 
is  important  to  note  that  the  remainder  of  Zalesak’s  procedure  must  still  be  carried  out. 
Antidiffusive  fluxes  acting  in  concert  on  cells  that  are  not  directional  extrema  of  p  could 
conspire  to  produce  new  local  maxima  or  minima  in  pn . 

The  resulting  limiter,  here  expressed  for  a  two-dimensional  problem  for  definiteness, 
consists  of  the  following  steps: 

A.  Establish  allowed  extrema  pmin  and  pmax  in  each  cell, 

ptf  =  mi*1 

Ptr  =  max  (pi,j-l’  Pi-hj'  Pij’  Pi+l,j’Pi,j+l)  ■ 

B.  Prelimit  the  antidiffusive  fluxes  along  each  coordinate  direction,  to  prevent  the  cre¬ 
ation  and  enhancement  of  directional  extrema  and  to  reverse  any  antidiffusive  flux 
that  is  directed  down  the  gradient  of  pe, 


F-+ 1/2,3  =  °i+W  max 


0,  min 


T?a 


icri+l/2,j  ViJ  Api_  1/2, 

&i+l/2,j  Vi+l,j  APi+3/2,j^J 


(11) 


8 


F?,j+i/2  =  °i,j+ 1/2  max 


0,  min^ 


rpa 

*ij+ 1/2 


Fi,j +1/2  i/2, 

ai,j+l/2  VSj+i  ^PiJ+3/2 


C.  Compute  the  total  incoming  and  outgoing  antidiffusive  fluxes  Fm  and  jF0Ut  i 
cell, 

—  max  (*7-1/2, j>°)  -  min  (*7+1/2,j>°) 

+  max  (f<_1/2,  o)  -  min  (i^'+1/2,  o)  , 

/  \  /  \ 


(12) 


in  each 


F°f  —  max 


T?a'  (Vl  _  min  ( TPa\  .  I 


(13) 


A  ^  .+1/2,*  >  i— l/2,j’  y 

+  max  (f<+i/2,  o)  -  min  o)  . 

D.  Determine  the  fractions  of  the  incoming  and  outgoing  fluxes  tha 
each  cell, 


\  / 

of  the  incoming  and  outgoing  fluxes  that  can  be  applied  to 

fij  =  (ptr  -  &)  v‘.i/F& 

f°T  =  (pl,i  -  <jn)  Vij/Ftf. 

).  Limit  each  antidiffusive  flux  so  that  it  creates  neither  an  undershoot  in  the  cell  it  i 
leaving  nor  an  overshoot  in  the  cell  it  is  entering, 


(14) 


is 


T7IC  _  77M 


Fc . 


*  x(  min  (fff,  flUj,  1),  H 

i+1/2j  \min(y53,/SSj,l),  o 

(  — (  -Tout  ^in  i  \  u 


if  *7+1/2,*  >  0; 

otherwise; 


*  i  i  /  r» 


This  procedure  differs  from  Zalesak’s  original  only  in  the  crucially  important  prelimiting 
step  B. 


6.  EXAMPLES 

Two  hydrodynamical  problems  will  be  used  to  illustrate  the  performance  of  the  pro¬ 
posed  multidimensional  FCT  flux  limiter.  The  first  is  a  two-dimensional,  Cartesian  muzzle- 
flash  simulation.  A  3.2  cm  x  3.2  cm  section  of  the  tube  is  represented  by  a  320x320  grid 
of  uniformly  space  cells.  Vertically  bisecting  and  extending  1.2  cm  into  the  domain  from 
the  left  edge  is  a  solid  divider  whose  thickness  is  .16  cm  (120x16  cells).  A  stoichiometric 
mixture  of  hydrogen  and  oxygen  gas,  whose  molecular  weight  is  12  and  ratio  of  specific 
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Fig.  3.  Results  of  a  Cartesian  muzzle- flash  simulation  are  shown  at  a  fixed  time  after  a  diaphragm, 
extending  from  the  end  of  the  horizontal  divider  to  the  bottom  of  the  domain,  bursts.  The  mass  density  is 
shaded  over  its  full  range  of  variation  (a  factor  of  5.5).  Zalesak’s  flux  limiter  was  used.  Note  the  prominent 
ripples  at  the  edges  of  the  light-  and  dark-shaded  regions,  behind  the  shock  front. 

heats  is  7/5,  fills  the  tube.  Behind  a  diaphragm  extending  from  the  end  of  the  solid  wall 
down  to  the  bottom  edge  of  the  domain,  the  mixture  is  at  a  pressure  of  20  atm  and  a  tem¬ 
perature  of  1000  K.  Elsewhere  the  pressure  and  temperature  are  at  standard  conditions. 
The  diaphragm  bursts  at  time  t  =  0,  allowing  the  hot,  pressurized  gas  to  flow  into  the 
chamber. 

The  Euler  equations  were  solved  using  a  two-dimensional  FCT  hydrodynamics  scheme, 
LCPFCT2,  described  elsewhere  (DeVore  1989).  It  is  the  2D  extension  of  Boris  and  Book’s 
one-dimensional,  low  phase  error  algorithm  LCPFCT  (Boris  &  Book  1976b;  Boris  et  al. 
1993),  and  shares  with  it  fourth-order  accuracy  in  phase  and  amplitude  at  long  wavelengths. 


*56,  N't  ^TSMSi.  .  3 


Fig.  4.  Same  as  Fig.  3,  but  using  the  new  monotone  flux  limiter.  The  prominent  ripples  evident  in  Fig. 
3  are  largely  eliminated. 

The  boundaries  were  closed  with  reflecting  conditions  at  the  top  and  bottom  of  the  domain 
and  along  the  tube’s  divider,  and  open  with  zero-gradient  conditions  at  the  left  and  right 
edges  of  the  domain.  A  Courant  number  of  .25  was  used. 

Computed  solutions  for  the  mass  density  after  an  elapsed  time  of  15  /xs,  using  Zalesak’s 
flux  limiter  and  the  proposed  monotone  limiter,  are  shown  in  Figs.  3  and  4,  respectively. 
Along  the  bottom  edge  of  the  domain,  the  solution  is  approximately  that  of  the  one¬ 
dimensional  Riemann  problem.  The  outermost  discontinuity  in  the  density  is  the  shock 
wave,  propagating  toward  the  right  edge  of  the  domain  near  the  bottom.  It  is  followed 
by  the  contact  discontinuity,  which  is  attached  to  the  open  end  of  the  divider  and  reaches 
about  75%  of  the  way  across  the  domain  at  the  bottom.  A  rarefaction  wave  propagates 
toward  the  left  below  the  divider.  The  expansion  of  the  fluid  upward  and  to  the  left  at  the 


Fig.  5.  Second  derivatives  along  the  horizontal  coordinate  of  the  mass  densities  obtained  with  Zalesak  s 


limiter  and  displayed  in  Fig.  3  are  shown.  The  color  table’s  range  is  restricted  to  saturate  the  shock  and 
contact  surface  and  bring  out  the  lower-amplitude,  numerical  fine  structure.  Note  the  rapid  oscillations 
between  maxima  and  minima  behind  both  the  shock  front  and  the  contact  surface. 

open  end  of  the  divider  creates  a  vortex  there.  These  principal  features  of  the  calculation 
are  well  defined  in  both  cases.  Zalesak’s  limiter,  however,  also  produces  a  substantial 
amount  of  oscillatory  fine  structure  in  its  solution.  This  is  evident  as  fingers  emanating 
from  the  shaded  bands,  especially  in  the  region  between  the  shock  wave  and  the  contact 
surface.  This  fine  structure  is  predominantly  numerical  in  origin. 

In  order  to  bring  out  the  spurious  structure  and  highlight  the  differences  between  the 
two  calculations,  the  second  derivative  of  the  mass  density  along  the  horizontal  coordinate 
is  displayed  in  Figs.  5  and  6  for  Zalesak’s  limiter  and  the  new  limiter,  respectively.  A  low 
value  for  the  extremal  amplitudes  spanned  by  the  color  table  was  selected,  so  that  the  shock 
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Fig.  6.  Same  as  Fig.  5,  but  obtained  from  the  density  field  in  Fig.  4,  which  used  the  monotone  flux  limiter. 
The  rapid  oscillations  clearly  evident  in  Fig.  5  are  greatly  suppressed. 

wave,  contact  surface,  and  edge  of  the  divider  strongly  saturate  the  figures,  allowing  the 
smaller-amplitude  fine  structure  to  be  brought  out.  Prominent  striations  fill  much  of  the 
volume  between  the  shock  and  contact  surface  and  to  a  lesser  extent  the  volume  inside  the 
contact  surface,  in  the  calculation  using  Zalesak’s  limiter.  These  features  are  much  weaker 
through  most  of  the  domain  when  the  new  limiter  is  used,  as  can  be  seen  by  comparing 
Figs.  5  and  6.  An  exception  is  the  small-amplitude  density  jump  between  the  shock  front 
and  the  contact  surface  that  is  evident  in  the  right  panel.  This  structure  is  faithfully 
reproduced  in  an  operator-split  integration  of  the  Euler  equations  for  this  problem  using 
LCPFCT  (not  shown).  It  apparently  is  an  artifact  of  the  initial  conditions,  developing  in 
the  early  evolution  of  the  system  as  the  shock,  contact  surface,  and  rarefaction  wave  form 
and  separate  from  their  common  origin  (Boris  &  Oran  1995).  Multiple  such  structures  can 
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be  seen  in  Fig.  5,  near  the  bottom  edge  of  the  domain  between  the  shock  and  the  contact 
surface.  These  additional  waves  evidently  are  excited  by  the  action  of  the  non-monotone 
flux  limiter. 

The  second  example  is  a  three-dimensional  simulation  of  a  subsonic  rectangular  jet 
(Grinstein  1995).  An  initially  laminar  air  jet  at  standard  temperature  and  pressure  issues 
from  a  rectangular  nozzle  of  aspect  ratio  4:1  at  Mach  0.6.  The  quiescent  background  gas 
also  is  at  STP.  In  order  to  focus  on  the  dynamics  of  the  individual  vortex  rings  of  the  jet, 
an  isolated  ring  is  puffed  out  by  closing  off  the  inflow  boundary  after  a  finite  time  equal 
to  the  ratio  of  the  jet  equivalent  diameter  to  the  flow  speed.  The  evolution  of  the  ring  is 
then  followed  as  it  convects  downstream. 

The  numerical  model  uses  LCPFCT2  to  solve  for  the  flows  in  the  2D  cross-stream 
plane,  and  LCPFCT  for  the  streamwise  direction.  A  150x110x110  grid  was  used  to 
represent  the  domain,  with  the  cells  evenly  spaced  in  the  shear-flow  region  of  the  jet 
and  geometrically  stretched  in  the  cross-stream  directions  outside.  The  Courant  number 
was  0.5.  Fixed  mass  density  and  velocity  conditions  were  specified  on  the  boundary  cells 
defining  the  jet  orifice,  with  free-slip  conditions  obtaining  elsewhere  on  the  entrance  plane. 
Advection  conditions  were  imposed  on  those  quantities  at  the  outflow  boundary,  and  all 
variables  satisfy  stagnation-flow  conditions  at  the  cross-stream  boundaries.  The  pressure 
satisfies  the  inviscid  ID  pressure  equation  at  the  jet  orifice  and  a  non-reflecting  condition 
at  the  outflow  boundary. 

Time  sequences  of  isosurfaces  of  the  vorticity  magnitude  are  shown  in  Figs.  7  and  8, 
obtained  using  Zalesak’s  limiter  and  the  monotone  limiter,  respectively,  in  the  2D  FCT 
module.  In  the  figures,  time  increases  from  left  to  right  and  from  bottom  to  top.  The 
bottom-most  six  frames  in  each  figure  show  the  self-induced  deformation  and  axis  switching 
of  the  vortex  ring.  The  highly  curved  corners  accelerate  ahead  of  the  ring  sides  and  toward 
the  centerline,  pulling  the  minor  axis  sides  along  with  them.  This  process  bends  the 
ring  along  its  major  axis.  The  increasing  curvature  at  the  midpoints  of  the  major  sides 
accelerates  those  portions  streamwise  toward  the  leading  minor  sides  but  away  from  the 
jet  centerline.  This  results  in  a  nearly  planar,  axis-switched  configuration  of  the  vortex 
ring  at  frame  6.  This  early  evolution  is  essentially  identical  in  the  two  cases,  save  for  some 
intermittent,  small-scale,  numerical  features  evident  with  Zalesak’s  limiter. 

Subsequently,  the  ring’s  new  major  sides  pinch  together  and  reconnect,  forming  a 
pair  of  vortex  rings  linked  on  the  underside  by  two  thin  threads.  This  is  shown  in  the 
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Fig.  7.  Self-deformation,  reconnection,  and  subsequent  merging  of  an  isolated  vortex  ring  are  shown.  The 
ring  is  puffed  from  a  Mach  0.6  air  jet  whose  aspect  ratio  is  4:1.  Isosurfaces  of  the  vorticity  magnitude  are 
shown  at  15%  of  the  peak  initial  vorticity.  Time  increases  from  left  to  right  and  from  bottom  to  top  in 
the  figure.  These  results  were  obtained  using  Zalesak’s  flux  limiter. 

following  six  panels  of  Figs.  7  and  8.  While  both  simulations  clearly  show  the  bifurcation 
of  the  ring,  the  fine  structure  associated  with  the  threads  bridging  the  two  daughter  rings 
increasingly  differs  between  the  two  limiters.  Zalesak’s  limiter  permits  fluctuations  on 
the  threads  of  vorticity,  which  show  up  as  spikes  attached  to  the  isosurfaces  and  lead  at 
frame  12  to  the  fragmenting  of  the  threads.  Note  that  the  two  daughter  vortex  rings 
appear  to  be  irrevocably  separating.  In  contrast,  with  the  new  limiter  the  fine  structure 
is  cleanly  and  clearly  represented  through  frame  12,  the  vorticity  threads  remain  intact, 
and  the  daughter  rings  stay  in  close  proximity  to  one  another.  Indeed,  as  the  assembly 
proceeds  further  downstream,  the  vortex  rings’  outer  edges  accelerate  ahead  while  their 
inner  edges  migrate  towards  the  centerline,  merge,  and  reconnect.  This  is  shown  in  the  top 
four  panels  of  Fig.  8,  and  in  even  later  images  in  Grinstein  (1995).  Features  at  still  smaller 


Fig.  8.  Vortex  ring  evolution  depicted  in  Fig.  7,  but  here  calculated  with  the  montone  flux  limiter.  The 
lower  twelve  frames  coincide  precisely  in  time  with  the  frames  of  the  preceding  figure. 

scales  than  before,  within  the  vortex  rings  themselves,  become  visible  and  can  be  tracked 
readily  during  this  phase  of  the  evolution.  This  example  dramatically  demonstrates  that 
the  dynamical  behavior  of  the  system,  above  and  beyond  the  aesthetics  of  the  simulation 
results,  can  be  influenced  significantly  by  the  choice  of  flux  limiter. 


7.  DISCUSSION 


An  improved  prescription  for  limiting  the  fluxes  in  multidimensional,  flux-corrected 
transport  calculations  has  been  proposed.  This  modification  to  Zalesak’s  (1979)  original 
multidimensional  limiter  is  motivated  by  the  observation  of  numerical  ripples  of  modest 
amplitude  in  two-  and  three-dimensional  FCT  solutions  of  hydrodynamical  and  MHD  sys¬ 
tems.  The  occurrence  of  such  ripples  is  traced  to  the  absence  of  a  monotonicity  constraint 
in  Zalesak’s  limiter,  which  nevertheless  is  positivity  preserving  and  prohibits  the  creation 
of  new  and  the  enhancement  of  existing  extrema.  The  recommended  modification  is  the 
addition  of  a  prelimiting  step  which  imposes  the  same  prohibition  on  directional  extrema, 
i.e.,  any  maxima  and  minima  of  the  transported  quantity  when  considered  along  each 
coordinate  direction  independently.  This  is  accomplished  by  applying  Boris  and  Book’s 
(1973)  original,  one-dimensional  flux  limiter  to  the  unidirectional  antidiffusive  fluxes. 

Two  numerical  examples  shown  illustrate  the  prevalence  of  small-scale,  relatively  low- 
amplitude  fluctuations  in  solutions  obtained  with  Zalesak’s  limiter,  due  to  the  lack  of 
a  monotonicity  constraint.  These  fluctuations  are  very  effectively  suppressed  by  the  new 
limiter,  leading  to  obvious  aesthetic  improvements  in  the  FCT  solutions.  More  importantly, 
the  rectangular  jet  simulation  demonstrates  compellingly  that  qualitative  changes  in  the 
system’s  dynamical  evolution  can  ensue  from  the  improved  representation  of  its  small-scale 
structures. 

Some  final  remarks  are  offered  concerning  one  other  facet  of  Zalesak’s  FCT  limiter: 
the  flexibility  available  in  establishing  the  allowed  extrema  in  each  computational  cell  (step 
A  of  the  limiting  procedure).  In  order  to  counter  FCT’s  tendency  to  ‘clip’  the  maxima 
and  minima  of  a  transported  quantity  and  form  three-cell- wide  extrema,  Zalesak  suggested 
looking  back  to  the  starting  values  at  each  timestep,  or  even  linearly  extrapolating  to  find 
extrema  between  cells,  in  addition  to  considering  the  positive  and  monotone  provisional 
solution  to  which  the  corrected  antidiffusive  fluxes  are  added.  This  generalization  seems 
to  be  advantageous  in  linear  advection  problems,  where  the  velocity  is  fully  specified  and 
minimizing  the  residual  numerical  diffusion  is  desirable.  In  solutions  of  nonlinear  systems 
such  as  the  Euler  equations  of  hydrodynamics,  on  the  other  hand,  this  increased  flexibility 
in  the  limiter  tends  to  create  more  numerical  ripples  of  greater  amplitude.  Fluctuations 
in  the  pressure,  e.g.,  create  associated  fluctuations  in  the  velocities,  and  the  resulting 
feedback  loops  are  detrimental  to  the  smoothness  of  the  solutions  obtained.  Restricting 
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the  determination  of  the  allowed  extrema  to  the  neighboring  provisional  values  along  the 
coordinate  axes  provides  sufficiently  accurate,  but  also  monotone,  solutions. 
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